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s  ABSTRACT 

The  fluid  mechanical  processes  which  characterize  a  transition  from 
deflagration  to  detonation  in  granular  beds  of  solid  propellant  have  not  at 
present  been  sufficiently  refined  to  allow  accurate  modelling  of  the  pheno¬ 
menon.  In  an  attempt  to  improve  this  situation,  this  report  has  investiga¬ 
ted  what  might  be  considered  basic  mechanisms  and  consequences  that  arise 
from  a  set  of  assumptions  for  the  governing  and  constitutive  equations  which 
take  into  account  the  two  phase  nature  of  this  problem.  A  qualitative  de¬ 
scription  of  the  flow  process  is  made,  based  on  observations  obtained  from 
DDT  experiments.  From  this,  certain  conclusions  are  reached  as  to  the  pro¬ 
perties  needed  by  propellants  to  exhibit  a  deflagration  to  detonation  transi¬ 
tion  (DDT).  The  numerical  integration  scheme  itself  is  examined  in  detail 
in  order  to  further  understand  the  consequences  of  its  use.  Also,  two 
scenarios  for  DDT  are  presented  which  exhibit  characteristics  similar  to 
those  derived  from  experimental  evidence.  Conclusions  as  to  the  direction 
of  future  research  are  made  based  on  the  results  obtained  from  the  work 
which  led  to  these  two  basic  mechanisms  for  DDT. 
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CHAPTER  ONE 


DESCRIPTION  OF  THE  FLUID  DYNAMICS  OF  DDT 

1 . l  Introduction 

The  work  outlined  in  this  report  is  a  continuation  of  the  research 
being  conducted  under  the  direction  of  Dr.  Herman  Krier  to  determine  a 
mechanism  which  explains  deflagration  to  detonation  transition  (DDT)  in 
porous  reactive  solids.  The  specific  regime  under  consideration  is  the  DDT 
phenomenon  observed  in  tightly  packed  beds  of  finely  granulated,  highly 
energetic  solid  propellants,  a  description  of  which  can  be  found  in  Krier 
and  Gokhale  [1]  with  further  results  to  be  found  in  Krier  and  Kezerle  [2] 
and  Krier,  Gokhale,  and  Hoffman  [3], 

The  method  used  in  analyzing  this  problem  is  basically  the  same  as 
that  developed  in  Krier  and  Kezerle  [2],  namely  the  use  of  six  conserva¬ 
tion  equations  based  on  the  concept  of  separated  continuum  flow  in  order  to 
accurately  solve  the  unsteady  fluid  dynamics  in  a  packed  bed.  The  exact 
form  of  the  equations  used  will  be  presented  later  in  this  report.  Previous 
to  the  use  of  this  approach,  a  system  of  equations  which  could  be  used  to 
analyze  the  same  problem  was  developed  by  Krier  and  Gokhale,  based  on  the 
continuum-mixture  theory  put  forth  by  Soo  [4].  Under  this  assumption,  the 
momentum  and  energy  equations  of  the  mixture  described  a  continuum  fluid  in 
terms  of  the  mean  mass  weighted  mixture  velocity  and  the  mixture  energy. 
However,  the  assumption  that  these  mixture  equations  describe  a  continuum 
(a  Navier-Stokes  fluid)  cannot  always  be  adequately  justified.  Thus,  sub¬ 
sequent  work  by  Krier  and  Kezerle  modified  this  approach,  and  assumed  that 

! 
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each  phase  independently  made  up  a  continuum.  The  important  difference 
that  can  now  be  seen  in  the  two  sets  of  equations  is  the  appearance  of 
inertial-coupling  (diffusion- like)  terms  in  the  balanced  equations  of 
momentum  and  energy.  (See  Soo  [5].) 

In  addition  to  these  conservation  equations,  state  equations  for 
both  the  gas  phase  and  the  solid  phase  of  the  bed  will  be  used  to  make  the 
computer  simulation  as  realistic  as  possible.  In  this  case,  the  gas  state 
equation  has  been  formulated  to  provide  accurate  results  in  the  nonideal, 
high-pressure  flow  which  is  under  consideration.  The  use  of  a  state  equa¬ 
tion  for  the  solid  phase,  which  represents  an  improvement  over  the  incom¬ 
pressible  assumption  used  in  previous  works  will  allow  variations  in  the 
particle  density  to  take  place. 

Due  to  the  separated  nature  of  the  flow,  there  must  be  supplemental 
constitutive  laws  which  account  for  the  interaction  of  the  two  phases  of 
the  flow.  These  laws  take  into  account  such  factors  as  an  interphase  heat 
transfer  and  a  gas-particle  drag.  In  the  studies  reported,  Krier-Gokhale  [1] 
and  Krier-Kezerle  [2],  a  compaction  resistance  law  was  used  that  later 
proved  to  be  unsatisfactory  when  subsequent  calculations  were  made  for 
packed  beds  of  increasing  length  [3).  This  particular  problem  was  solved  by 
linking  the  compaction  resistance  (particle  stress)  to  the  amount  of  force 
being  applied  to  the  bed  by  the  interphase  drag.  Each  individual  phase  also 
has  a  set  of  constitutive  relations  which  are  necessary  to  completely  de¬ 
scribe  the  flow.  These  include  such  equations  as  a  law  describing  the  axial 
particle  stress  (related  to  the  compaction  resistance  law),  a  model  to  de¬ 
scribe  the  movement  of  this  stress  wave,  a  pressure-dependent  burning  rate 
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law  for  the  particles,  and  a  temperature-dependent  gas  viscosity  law. 
Variations  in  any  one  of  these  constitutive  laws  caused  its  own  change 
in  the  progression  of  the  flow  through  the  bed.  Again,  results  from 
variations  in  some  of  these  laws  will  be  presented  later. 

1 . 2  bef lagration-to-Detonation  Transition  (DDT) 

As  discussed  in  Krier  and  Kezerle  [2],  experimental  studies  of  the 
spontaneous  deflagration-to-detonation  transition  (DDT)  in  a  packed  bed  of 
granular  or  solid  propellant  indicates  that  in  the  steady-state  detonation 
phase,  a  pressure  shock  wave  precedes  the  flame  (or  ignition)  front.  In 
the  build-up  to  this  phase,  there  appears  to  be  two  scenarios  which  describe 
the  pre-detonation  fluid  dynamics,  both  of  which  can  be  inferred  from  ex¬ 
perimental  evidence.  The  first  mechanism  results  from  the  work  by  Ber- 
necker  and  Price  [6] ,  which  assumes  that  the  pressure  front  (not  yet  a 
shock  wave)  is  initially  behind  the  flame  front.  But  due  to  the  rapid 
generation  of  gases  by  the  burning  propellant,  the  pressure  front  is 
simultaneously  increased  in  magnitude,  steepened  in  its  gradient,  and 
accelerated  towards  the  flame  front.  The  point  at  which  the  pressure 
front  overtakes  the  flame  front  is  then  defined  as  the  transition  point, 
after  which  the  flame  front,  now  preceded  by  the  pressure  front,  moves 
down  the  bed  in  a  self-sustaining  detonation  wave  at  a  speed  greater  than 
the  deflagration  wave.  (See  Fig.  1.1.) 

The  second  DDT  mechanism  can  be  traced  to  work  first  done  by  Matek 
[7]  for  DDT  in  solid  explosives  which  was  later  expanded  upon  by  Tarver 
et  al ■  [8].  In  this  model,  the  rapid  build-up  of  gas  pressure  behind 
the  flame  front  sends  compression  waves  into  the  compacted,  but  as  yet 


experiments  of  Bernecker  and  Price  [6]. 


unignited  propellant.  Due  to  the  increase  in  sound  speed  behind  suc¬ 
cessive  compression  waves,  each  wave  will  eventually  overtake  its  pre¬ 
decessor,  which  will  result  in  the  formation  of  an  insipient  shock 
dig.  i.2j.  The  shock  is  assumed  to  form  some  distance  ahead  of  the 
flame  front,  resulting  in  a  steady-state  detonation  wave  moving  forward 
through  the  unignited  portion  of  the  bed,  as  well  as  a  retonation  wave 
moving  back  towards,  and  eventually  meeting,  the  flame  front.  This  model 
has  been  given  some  experimental  support  for  DDT  in  granular  beds  from  work 
by  Calzia  and  Carabin  (  9].  (See  Fig.  1.3.)  It  should  be  noted  at  this 
point  that  the  first  mechanism  could  be  considered  a  special  case  of  the 
second  if  both  the  flame  front  and  gas  pressure  front  were  to  arrive  at 
the  point  of  formation  of  the  insipient  shock  at  the  exact  time  of  shock 
format  ion. 

1 . 3  Areas  of  Invest  _i£at  ion 

The  work  presented  in  Krier  and  Kezerle  [2]  developed  the  necessary 
conservation  equations  and  constitutive  laws  to  perform  the  pre-detonation 
fluid  dynamics,  which  is  a  necessary  part  of  either  of  the  previously 
mentioned  DDT  models.  Using  this  set  of  laws,  an  effort  to  find  the  DDT 
mechanism  based  on  the  first  mechanism  was  made  by  utilizing  a  selective 
cariation  of  parameters  in  this  set  of  equations.  The  results  of  this 
effort  seemed  to  indicate  that  unless  some  alternative  jump- like  function 
were  introduced  in  the  equations,  one  does  not  predict  the  formation  of  a 
reaction  front  at  speeds  in  the  range  of  10  mm/psec  and  being  supported  in 
advance  by  the  detonation  shock.  Some  candidates  which  include  the  physics 
to  allow  for  such  jump-like  function  are  the  following  concepts: 


Distance 
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Fig.  1.3  Detonation  to  deflagration  transition  mechanism  supported  by 
experiments  of  Caizia  and  Carabin  [91 . 


(1)  At  some  critical  high  pressure,  P* ,  or  at  some  critical  (dP/dt)*,  the 
material  burns  at  a  "super"  rate  different  from  the  rates  assumed  and 
extrapolated  from  the  lower  pressure  data  base. 

(21  Since  the  burning  rate,  r ,  appears  only  in  the  mass  generation  term, 

F,  it  should  logically  follow  that  one  should  look  at  other  terms  in 
F  for  use  as  a  discontinuous  type  function.  For  example,  since 
T  =  p^rfl  -  0) (3/tp) ,  one  might  wish  to  formulate  a  state  equation  in 
which  either  density,  ,  or  solids  loading,  (1  -  0) ,  increases  dis- 
continuously ,  as  some  high  (critical)  pressure. 

(3)  Again  looking  at  the  T  term,  one  could  assume  that  at  a  given  instant, 
a  critical  pressure  could  fracture  the  particles  thereby  drastically 

increasing  the  surface  to  volume  ratio,  (3/r  ). 

P 

txperimental  studies  to  investigate  these  possibilities  should  be 
investigated  even  though  verifying  the  results  would  probably  prove  very 
difficult.  A  limited  number  of  these  potential  jump-like  functions  have 
been  studied  analytically  and  are  reported  herein. 

With  regards  to  the  second  model  for  DDT,  the  codes  developed  and 
used  by  Krier-Cokhale  [1]  and  Kr i er- Kezerle  [2]  made  use  of  a  compaction 
resistance  law  that  was  a  function  of  only  the  porosity,  <J> ,  which  could  be 
interpreted  as  the  particle  stress.  Due  to  the  nature  of  the  equations, 
the  minimum  porosity  always  occurred  in  association  with  the  gas  pressure 
front  and  thus  the  stress  wave  would  not  propagate  into  the  bed  in  a  man¬ 
ner  similar  to  a  compression  wave,  but  merely  remain  tied  to  the  gas  pres¬ 
sure  front.  In  order  to  obtain  the  presumed  stress  wave  motion,  it  was 
determined  that  the  particle  stress  would  need  to  be  predicted  in  such  a 
way  as  to  correct  the  problem  encountered  earlier  in  the  longer  length 


beds  and  to  decouple  the  motion  of  the  stress  wave  from  the  motion  of  the 
gas  pressure  front.  This  was  accomplished  by  locating  the  point  of  maxi¬ 
mum  stress  (presumed  to  be  a  function  of  both  gas  pressure  and  the  drag 
force]  and  allowing  that  stress  to  propagate  forward  at  a  sound  speed 
associated  with  that  stress.  The  exact  form  of  this  sound  speed  relation 
will  be  detailed  in  Appendix  A.  Since  there  are  no  distinct  viscosity 
type  terms  associated  with  the  solid-phase  equations  which  would  normally 
account  for  the  motion  of  this  type  of  disturbance,  it  was  felt  that  im¬ 
posing  this  kind  of  assumption  on  the  stress  wave  would  simply  be  a  short¬ 
hand  way  of  describing  the  dynamic  compression  of  a  "homogeneous  solid". 

In  general  the  method  of  investigation  presented  here  is  constructed 
in  such  a  manner  so  as  to  provide  conditions  in  which  the  DDT  phenomenon 
can  manifest  itself.  This  is  done  by  assuming  that  a  fraction  of  the  pro¬ 
pellant  grains  are  ignited  at  one  end  of  a  closed  chamber.  The  mass 
generated  in  the  ignition  region  accelerates  the  resulting  hot  gases  forward, 
with  the  region  behind  these  hot  gases  rapidly  increasing  in  pressure.  This 
pressure  rise  is  due  to  the  large  quantities  of  gaseous  material  generated 
when  the  propellant  is  assumed  to  obey  a  pressure  sensitive  burning  rate 
law.  At  this  point  it  is  presumed  that  one  of  the  two  transition  mechanisms 
described  above  will  cause  the  bed  to  exhibit  a  DDT.  Regardless  of  the 
mechanism  which  causes  the  transition,  a  detonation  can  be  said  to  have 
begun  when  the  pressure  front  precedes  the  flame  front  and  both  move  through 
the  bed  at  a  speed  which  is  characteristic  of  a  detonation  wave  for  the  type 
of  propellant  and  initial  porosity  used.  In  this  report,  the  flame  front 
will  be  defined  as  the  locus  of  points  which  mark  the  position  in  the  bed 
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of  the  initial  particle  ignition.  The  definition  of  the  pressure  front 
will  depend  on  the  transition  mechanism  being  considered  and  thus  will  be 
deferred  until  Chapter  Three  in  which  particular  cases  will  be  discussed 
in  detail . 


CHAPTER  TWO 


THE  GOVERNING  EQUATIONS  FOR  UNSTEADY, 
ONE -DIMENSIONAL,  TWO-PHASE  FLOW 


2. 1  Introduction 

The  set  of  conservation  equations  and  constitutive  laws  that  make 
up  the  governing  equations  for  this  investigation  are,  for  the  most  part, 
the  same  as  those  found  in  Krier  and  Kezerle  [2],  As  mentioned  earlier, 
the  approach  taken  in  developing  the  conservation  equations  assumes  that 
there  are  two  distinct  continua,  one  for  solids  and  one  for  the  gas,  each 
moving  through  its  own  control  volume.  Due  to  this  approach  the  sum  of 
these  two  volumes  must  represent  an  average  mixture  volume  while  at  the 
same  time  the  equations  which  describe  the  two  continua  must  account  for 
the  effect  that  one  flow  has  on  the  other.  To  obtain  this,  distinct  equa¬ 
tions  for  continuity,  momentum,  and  energy  are  written  for  each  phase  which 
recognize  that  each  phase  occupies  only  part  of  the  total  volume  and  uti¬ 
lizes  inertial-coupling  terms  which  disappear  when  the  two  sets  of  equa¬ 
tions  are  summed  together.  A  detailed  derivation  of  these  equations  is 
presented  in  Appendix  A  of  Krier  and  Kezerle  [2].  However,  due  to  the  two 
detonation  mechanism  theories  being  considered,  certain  modifications 
must  be  included  and  will  be  noted  when  made. 

2. 2  Assumptions 

The  following  are  the  set  of  assumptions  and  limitations  placed  on 
the  governing  equations  used  in  this  analysis. 

(1)  Each  phase  of  the  flow  is  a  continuum  thus  allowing  for  unique  deriva- 
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(2) 


(3) 


(4) 

(5) 

(6) 

(7) 

(8) 

(9) 

(10) 

(11) 


(12) 


Although  each  phase  is  a  continuum,  they  are  considered  interdis- 
persed  which  is  reflected  in  the  conservation  equations  by  the 
presence  of  inertial-coupling  terms. 

The  flow  is  quasi-one-dimensional ,  that  is,  the  total  cross-sectional 
area  of  the  flow  is  equivalent  to  the  sum  of  the  cross-sectional  area 
of  each  phase. 

During  combustion,  the  solid  phase  always  loses  mass  to  the  gas 
phase  (i . e. ,  T  1  0) . 

The  equations  are  laminar  in  that  the  turbulence  resulting  from  the 
two-phase  nature  of  the  flow  has  been  averaged  out. 

All  gases  obey  a  nonideal  Nobel-Abel  equation  of  state  with  a 
variable  covolume. 

All  gases  are  inviscid  in  nature  with  the  exception  of  their  action 
in  the  drag  relation. 

In  those  relations  which  use  the  gas  viscosity  or  conductivity, 
these  properties  are  assumed  to  be  functions  of  temperature. 


The  specific  heats,  c  and  c  ,  of  both  phases  are  considered  constant. 

p  v 


The  solids  obey  the  Tait  equation  of  state. 
Conductive  or  radiative  heat  transfer  is  neglected. 
All  particles  are  spherical. 


2 . 3  Conservation  Equations 

In  order  to  uniquely  describe  the  properties  of  the  two-phase  flow, 

the  following  nine  variables  must  be  determined:  p,p,u,u,T,T, 

g  P  g  P  g  P 

Pg.  Pp»  and  <t.  Since  there  are  nine  unknowns,  nine  equations  must  be 
supplied  in  order  to  find  a  singular  solution.  The  following  conservation 
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equations  can  be  used  to  provide  six  of  the  necessary  relations. 
Gas  Continuity 

3tpiV  , 

3t  '  3x 


Solid  Continuity 

^2=  3(p2V  p 

3t  3x 


Gas  Momentum 


3(Plug) 

3t 


— tt — -  <(>  &  -  F  +  Tu 

dx  dx  p 


Solid  Momentum 


3(p2up) 

3t 


3 (p.u2)  3P 

-Ip"  -zf 


fu 


Gas  Energy 


’V1’ 

at 


3  (u  E  p .  +u  <J>P  ) 
g  g„  1  g^  e' 


g  g 


3x 


+  r 


p  CHEM 

.  2  g 


Fu 


Particle  Energy 


3(Ep  p2)  3(upEpTP2+up(l-<!,)Pp) 


at 


ax 


+r 


_E  +  echlm 

2  p 


+  Fu  +  Q 
P 


Where,  we  define 


Phase  densities 


Porosity 


P1  =  *Pg 


P-,  =  (i-4>)Pr 


v  v 

=  ;  solid  loading  :  C 1  -03  = 

*  u:  ..  ‘  w 


Mi  x 


Mix 


(2.1) 

(2.2) 

(2.3) 

(2.4) 

(2.5) 

(2.6) 

(2.7a) 

(2.7b) 


- 
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d. 


Total  internal  energy  (gas) 


E 


e 

g 


1  ? 
+  ■=■  uz 
2  g 


Total  internal  energy  (particles) 


E 

P 


and 


E 

g 


and 


E 

P 


c  T 
v  p 
P 


(2.7c)’ 


(2 . 7d) * 


In  this  case,  the  subscript  "g"  denotes  the  gas  phase  while  the  "p" 
indicates  the  particle  phase.  Here,  refers  to  the  chemical  energy 

released  in  burning.  The  seventh  and  eighth  equations  needed  are  the 
equations  of  state  associated  with  each  phase.  For  the  gas,  the  Nobel- 
Able  equation  is  used 


p  R  T 


1 g  1 -p  B 
g  v 


(2.8) 


where  B v  is  the  covolume,  a  term  that  is  a  function  of  gas  density.  The 
use  of  results  in  a  useful  nonideal  equation  at  high  gas  pressures.  In 
this  study  it  is  assumed  that 


1-p  B 
g  v 


a  +  bp  + 


cpz  +  dp3 
g  g 


(2.9) 


where  a  =  1.0,  b  =  1.0,  c  =  0.5,  d  =  0.3,  and  p  is  in  grams/cm3. 

g 

This  equation  results  from  an  assumption  that  the  gas  particles  behave  as 
hard  spheres  during  any  interaction.  Appendix  I)  will  show  that  one  can 
use  Eq.  (2.7c)  for  the  gas  internal  energy  given  the  form  of  Eqs.  (2.8) 
and  (2.9). 


Here  c  =  0.30316  BTTJ/Lbm°R  and  c  =  0.42442  BTU/Lbm"R.  Both  are 
Vg  VP 

assumed  average  constants,  independent  of  temperature. 
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2.4  Constitutive  Relation  for  P 
_ £ 

As  described  by  Wallis  [10],  a  packed  bed  placed  under  a  compressive 
load  can  be  further  compacted.  However,  there  is  a  force  which  resists 
this  compaction  that  depends  on  the  stress-strain  relationship  of  the  par¬ 
ticle  lattice  which  is  not  necessarily  the  same  as  that  of  a  homogeneous 
solid  made  from  the  same  material.  The  compressive  load  on  the  bed  will 
be  split  between  the  two  phases  in  proportion  to  the  porosity.  So  the 
resultant  force  on  the  particles  will  be  a  function  of  the  porosity,  the 
porosity  gradient,  and  possibly  other  factors.  There  can  be  a  variety  of 
formulations  that  could  be  used  to  relate  this  particle-particle  inter¬ 
action  through  a  stress,  which  is  termed  here  as  P  .  That  is, 

P 

P  =  P  (<(>,  (AV  )...  etc.) 
p  p  T  dx  Mix  1 

One  formulation  would  be  to  assume  that  P^  is  a  function  of  porosity 
only.  This  approach,  taken  by  Kuo  and  Summerfield  [11],  resulted  in  an 
equation  of  the  form 


B[(l-*c)_1  -  (l-'*’)'1] 

(T^S 


if 


if 


C2.10) 


where 


is  the  porosity  above  which  the  particles  do  not  touch.  (See 


Appendix  B.)  In  order  to  get  some  idea  of  the  validity  of  this  equation, 
the  relation  for  P^  can  be  used  to  derive  a  bulk  modulus  for  the  particle 
lattice  which  one  would  expect  to  increase  as  the  solids  loading  increases. 
One  defines  the  bulk  modulus,  K,  as 
3P 


K  = 


3VMix 


VMix 


(2.11) 
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which  for  Eq.  (2.10)  becomes 


K  = 


j-B  [(l-ctg-'ll-*)'2 


2(l-4>)-3] 


L±  1  V 
VMix  j  Mix 


(2.12) 


In  the  above  relation  it  can  be  shown,  that  for  a  negligible  change  in  the 


volume  of  the  solids,  that  d<f>/dVMix  =  (1  -4>) /vMix 

As  indicated  in  Fig.  2.1  the  bulk  modulus  does  not  increase  as 
porosity  decreases  and  is  therefore  suspect.  In  addition,  this  formulation 
does  not  take  into  account  certain  constraints  which  are  present  if  one 
considers  only  spherical  particles.  In  this  case  the  particle-particle 
interaction  term  must  prevent  compaction  below  a  porosity  of  (J)  =  0.2595. 
This  limitation  results  from  the  fact  that  the  greatest  possible  compaction 
of  spherical  particles  results  from  a  face  centered  cubic  lattice,  where 


^MIN  1 


0.2595 


(2.15) 


(Details  of  this  calculation  and  other  limits  on  the  porosity  of  sphere  can 
be  found  in  Appendix  B.)  Given  this  type  of  constraint,  there  is  no 
guarantee  that  Eq.  (2.10)  would  under  all  conditions  provide  such  a  proper 
lower  packing  limit. 

One  approach  toward  satisfying  this  limit  is  to  modify  the  particle 
momentum  equation  so  that  the  particles  cannot  be  accelerated  forward  once 
critical  lower  porosity  level  is  reached.  A  modification  of  this  type 
would  cause  Eq.  (2.4)  to  look  like 


3(p2y  =  _  a(p2up) 

3t  3x 


—  ■ 

3P 

l  -  f(4>) 

F  -  (l-<fr) 

_ 

L_ 

ru  + 

p 


(2.14) 


Bulk  Modulus ,  K  (kpsi) 


Fig.  2.1  Bulk  modulus  versus  porosity  for  particle-particle  interaction 
law  of  Kuo  and  Summerfield  [11]. 
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where  f(<j>)  is  defined  as 


f  (<J>)  = 


1 


0 


♦ 


LL 


Z 


1.0 


4> 


>  *UL 

UL  "  ♦  -  *LL 
<  *LL 


(2.15) 


Here  <(>  is  the  upper  limit  on  the  porosity  for  a  packed  bed, 

lower  limit  on  the  porosity,  and  Z  is  an  exponent  chosen  to  reflect  the 

properties  of  the  particle  lattice.  Figure  2.2  shows  the  variation  of 

f((f>)  with  porosity  for  various  values  of  the  constant,  £. 

Equation  (2.14)  indicates  that  the  drag  force  has  a  decreasing 

effect  on  the  particle  acceleration  as  some  lower  limit  to  the  porosity 

is  reached.  However,  the  gas  phase  will  still  feel  the  full  effect  of 

the  drag  on  its  deceleration.  (Recall  that  most  experiments  used  to 

determine  the  drag  relation  for  a  porous  media  are  made  through  a  bed 

of  essentially  non-moving  particles.)  The  apparent  discrepancy  in  these 

two  statements  is  the  origin  of  the  constitutive  relation  for  P  .  This 

P 

relation  comes  about  through  the  assumption  that  the  drag  force  not  in¬ 
volved  in  accelerating  the  particles  is  instead  stressing  the  particle 

bed,  which  is  the  definition  for  P  .  Thus,  P  is  now  related  to  the 

P  P 

drag,  F,  in  the  following  manner,  namely 


P 

P 


(F) 


(2.16) 


An  important  observation  about  this  type  of  formulation,  i.e., 
either  Eq.  (2.10)  or  (2.16),  is  that  the  stress  on  the  particles  is 
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localized.  The  particles  are  not  dynamically  loaded  as  would  happen  if 
the  solid  were  cemented  into  a  pseudo-homogeneous  material.  In  this 
latter  case,  fundamental  continuum  mechanics  predicts  that  a  stress  at 
one  point  in  the  material  would  propagate  into  an  unstressed  region  at 
a  rate  which  is  proportional  to  both  the  average  sound  speed  and  the 
particle  velocity  (i.e.  ,  a  right  running  or  left  running  characteristic). 
Since  the  particles  are  now  considered  compressible,  the  sound  speed 
will  vary  proportionately  with  the  solid  phase  density,  the  details  of 
which  are  presented  in  Appendix  A. 

Therefore,  an  alternate  constitutive  relation  (if  one  wishes  to 
include  dynamic  compressibility)  is  to  model  the  particles  as  a  contacting, 
but  obviously  very  porous,  solid  and  to  assume  that  while  the  particles 
are  in  contact  with  each  other  a  stress  wave  whose  magnitude  is  defined  in 
Eq.  (2.16)  moves  through  the  bed  at  a  speed  equal  to  the  sum  of  the 
sound  speed  behind  the  wave  and  the  local  particle  velocity. 


2. 5  Additional  Constitutive  Relations 

In  order  to  complete  the  analysis  of  Eqs.  (2.1)  through  (2.6)  it  is 
necessary  to  define  certain  criteria  and  interaction  laws  with  the  use  of 
known  physical  relationships  or  through  the  use  of  experimentally  estab¬ 
lished  laws.  These  include 

(1)  An  ignition  criterion  based  on  the  bulk  temperature  of  the 
particles. 

(2)  The  propellant  burning  rate,  r  . 

(3)  The  gas  generation  rate,  T,  which  for  spherical  particles  is 


r  = 


lrp 


d-^Ppfp 


(2.17) 
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(4)  The  interphase  heat  transfer  rate,  Q,  which  is  defined  as 


Q  =  h  (T  -T  )  C4>)  (S/V) 
Pg  g  P  P 


(2.18) 


where  (S/V)p  is  the  surface  to  volume  ratio  for  the  particles 

h  Is  the  interphase  heat  transfer  coefficient,  which  is  dis- 
Pg  v 

cussed  in  Appendix  C. 

(5)  The  interphase  drag,  F,  which  is  defined  as 


F  = 


p(U  -U  )(1-<|>)2 


4r2 

P 


Pg 


(2.19) 


where  p  is  the  gas  phase  viscosity  and  f  is  the  interphase  drag 
coefficient  which  is  presented  in  Appendix  C. 

(6)  The  temperature  dependent  gas  viscosity. 

The  constitutive  relations  presented  here  differ  very  little  from  those 
presented  in  Krier  and  Kezerle  [2],  where  additional  background  has  been 
given. 


2.6  Numerical  Integration  Technique 

In  order  to  solve  this  set  of  hyperbolic  nonlinear  partial  differ¬ 
ential  equations,  the  two  step  predictor-corrector  numerical  integration 
scheme  developed  by  MacCormack  was  used.  This  is  the  same  type  of  scheme 
as  detailed  in  Ref.  2  with  a  few  exceptions  which  will  be  outlined  below. 

Throughout  this  investigation  the  total  bed  length  was  usually  set 
at  3  inches  with  the  Ax  set  at  0.05  inches.  This  value  of  Ax  was  de¬ 
termined  to  provide  the  largest  mesh  size  for  repeatable  calculations. 

The  value  of  L  =  3  inches  gave  moderate  computing  charges.  The  At 


(for  this  Ax)  is  calculated  using  the  Courant ,  Friedrichs,  Levy  stability 


criteria  [12],  namely 


At.Ll'  *.  J.u|]  <  i  (2.20) 

Ax 

Here,  C  and  |u|  are  the  maximum  gas  phase  or  solid  phase  sound  speed  and 
phase  velocity,  respectively,  and  not  the  mixture  average  values  as 
utilized  by  Krier  and  Kezerle  [2],  This  was  done  to  provide  more  stable 
solutions  for  all  classes  of  problems  studied  here.  To  satisfy  the  in¬ 
equality  in  Eq.  (2.20)  the  value  of  unity  was  set  equal  to  0.7. 

Even  with  the  use  of  a  guaranteed  stability  criterion  for  the 
linear  form  of  the  hyberbolic  equations,  there  were  conditions  under 
which  the  computer  code  would  become  unstable  for  our  nonlinear  set  of 
equations.  Although  it  was  found  that  these  instabilities,  caused  by 
severe  gradients  in  one  or  more  of  the  partial  differential  equations, 
could  be  removed  in  some  cases  by  utilizing  a  smaller  Ax  (and  its  associ¬ 
ated  reduced  At),  an  artificial  smoothing  routine  was  adopted  to  insure 
stability.  A  detailed  discussion  of  this  smoothing  technique  can  be 
found  in  the  Ph.D.  dissertation  of  S.  S.  Gokhale  [13].  The  addition  of 
this  smoothing  technique  did  not  take  place  until  later  in  this  investi¬ 
gation.  Thus,  most  of  the  results  presented  here  did  not  employ  arti¬ 
ficial  smoothing  and  hence  the  range  of  input  parameters  was  often  limited. 

Ihe  boundary  conditions  chosen  for  this  investigation  are  the  same 
as  those  used  by  Krier  and  Kezerle  with  the  exception  of  the  choice  of 
gradien*  at  the  wall.  In  this  work  the  code  was  programmed  to  insure 
that  all  gradient  terms  at  both  end  walls  are  zero  at  all  times. 
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Two  sets  of  initial  conditions  were  utilized  in  this  investigation, 
the  choice  of  which  one  to  use  being  somewhat  arbitrary.  The  first  set 
of  conditions  were  similar  in  form  to  those  used  by  Krier  ar.i  Keierle, 
namely. 


(1)  T  and  T  were  set  by  a  lineal  decay  over  the  first  five 

g  P 

grid  points. 

(2)  was  subject  to  an  exponential  decay  throughout  the  bed. 

(31  Velocities  of  both  phases  were  set  to  zero  over  the  entire 

length  of  the  bed. 

(4)  The  porosity  was  set  at  some  selected  value  which  was  uniform 
across  the  total  bed  length. 

The  second  set  of  initial  conditions  came  about  from  the  feeling 
that  it  was  inappropriate  to  model  a  situation  in  which  there  is  no  gas 
motion  in  the  presence  of  a  gas  pressure  gradient.  Thus,  the  second  set 
of  conditions  assumes  that  additionally  at  t  =0,  the  pressure  is  a 


uniformly  low  constant  value. 
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CHAPTER  THREE 

COMPUTED  RESULTS 


3 . 1  Introduct ion 

As  discussed  in  Chapter  One,  the  investigations  described  in  this 
report  consist  of  a  more  detailed  analysis  of  the  fluid  mechanical  proper¬ 
ties  leading  up  to  the  detonation  transition  of  a  porous  reactive  bed  as 
opposed  to  those  studies  carried  out  by  Kri er -Cokhale  [1]  and  Krier- 
Keterle  [2J.  Specifically,  changes  in  the  solid  phase  equations  used  by 
krier  and  Keterle  were  made  in  an  attempt  to  gain  a  better  understanding 
of  the  properties  of  a  bed  of  closely  packed  particles  and  to  determine 
what  conditions  are  required  for  this  bed  to  exhibit  a  DDT.  These  changes 
included  the  addition  of  a  solid  phase  state  equation,  the  definition  and 
use  of  a  particle  pressure,  and  modi  fications  to  some  of  the  various  con¬ 
stitutive  laws  which  are  necessary  to  completely  describe  the  flow.  (See 
Appendix  C.)  A  discussion  of  the  effects  of  each  of  these  changes  on  the 
unsteady  fluid  dynamics  in  the  bed  will  be  presented  below. 

However,  to  provide  a  comparison  with  the  calculations  presented  in 
the  earlier  works  of  kr ier-Cokha  le  and  Kr i er-Kecer le  ,  a  baseline  case  was 
constructed  to  closely  reflect  typical  results  found  in  those  studies, 
fable  3.1  lists  the  most  important  inputs  for  both  the  baseline  calculation 
and  for  the  subsequent  parameter  variation  calculations.  The  results  of 
these  calculations  will  be  presented  using  one  or  more  of  the  following 
graphical  formats:  pressure,  temperature  and  porosity  versus  distance,  and 
the  ignition  front  position  (flame  front!  versus  time. 
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TABLE  3.1 

TYPICAL  INPUT  DATA 


Parameter 

Value 

English 

SIU 

Initial  bed  temperature,  T  , 

T 

P 

530°R 

294°K 

Ignition  energy,  AEIGN 

165  BTU/LB 

m 

383  KJ/Kg 

Ignition  temperature,  T^^ 

545°R 

303°K 

Initial  bed  porosity 

0.4 

0.4 

Propellant  burning  rate 
constant,  b^ 

9.0xl0-4  in/sec-psin 

3. 3 2x1  O'7  cm/sec-Pa11 

Propellant  burning  rate 
index,  n 

0.90 

0.90 

Initial  propellant  density. 

PP 

0.0571  LB  /in3 
m 

1.581  gra/cm3 

Initial  bulk  modulus,  Kq 

2.00x10s  LBf/in 

1.38  GPa 

Constant  volume  specific  heat 
of  the  propellant,  cy 

VP 

0.30316  BTU/LB  °R 
m 

1.2665  J/gm°K 

Initial  grain  radius,  r^ 

4 . 0x10” 3  in 

0.1016  mm 

Chemical  energy  released, 

E  -  (E  -E  )CHEM 

CHEM  1  g  p-1 

2360.9  BTU/LB 

m 

5479.6  KJ/Kg 

Molecular  weight  of  the 
gas,  MW 

22.6  LB  /LB  . 

m  mole 

Covolume  of  propellant  gas. 

B 

V 

29.85  in3/LB 

m 

1.078  cm3/gm 

Specific  heat  ratio  of  gas. 

Y 

1.252 

1.252 

Constant  volume  specific 
heat  of  gas,  c 

Vg 

0.42442  BTU/LB  °R 
m 

1.7731  Joule/gm-°K 

Gas  viscosity,  ug 

2.49xl0‘6  LB  /in-sec 
m 

4.45xl0-4  gm/cm-sec 

Universal  gas  constant,  R 

1.9869  BTU/LB  .  °R 
mole 

8.3005  Joule/gm-°K 

Total  bed  length,  L0 

3.0  in 

76.2  mm 

Percent  heat  transfer  to 
particle  after  ignition 


10% 


10% 
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3. 2  Baseline  Case 

An  initial  baseline  case  was  constructed  that  closely  resembles 
the  conditions  studied  by  Krier  and  Kezerle  [2] .  This  calculation  assumes 
that  the  solid  phase  is  incompressible  and  that  the  packed  bed  obeys  the 
particle-particle  interaction  law  develed  by  Kuo-Summerfield  [11].  As  in 
Krier-Kezer le  an  initial  uniform  porosity  of  0.40  (solids  loading  of  60%) 
was  used. 

Upon  comparison  of  the  two  baseline  cases,  it  was  found  that  after 
a  certain  percentage  of  the  bed  had  been  ignited,  the  flow  properties  of 
each  were  roughly  the  same  in  magnitude.  However,  the  time  necessary  to 
reach  that  percentage  of  ignition  was  now  approximately  30  percent  larger 
when  compared  to  those  results  reported  in  Ref.  2.  This  was  due  mainly  to 
differences  in  the  value  specified  for  the  specific  heat  for  the  propel¬ 
lant  and  the  coefficient  in  the  interphase  drag  relation  currently  being 
used . 

The  development  of  the  gas  pressure  distribution  (Fig.  3.1)  obviously 
is  similar  to  that  presented  by  Krier  and  Kezerle  in  that  the  pressure 
front  shows  a  distinct  tendency  to  build  into  a  continental  divide.  A  con¬ 
tinental  divide  is  defined  here  as  a  gas  pressure  spike  which  usually 
appears  in  the  region  of  the  deflagration  front.  This  spike  is  significant 
in  that  the  pressure  gradients  which  have  now  been  formed  will  cause  the 
gas  to  move  both  towards  and  away  from  the  burning  zone.  The  reason  that 
such  a  divide  forms  is  due  to  the  continuing  availability  of  gas  volume  in 
the  section  initially  ignited  as  the  propellant  burns  away  in  an  ever  ac¬ 
celerating  manner  as  the  pressure  rises.  A  typical  value  for  this  gradient 


Fig.  3.1  Pressure  distribution  during  flame  spreading  in  an 
initially  packed  bed  of  small  particles  of  solid 
propellant.  (Figs.  3. 2-3. 4  show  variation  in  other 
flow  parameters  for  the  same  case  as  calculated  for 
Fig.  3.1.) 
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at  t  =  60  nsec  and  x  =  44.5  min  (1.75  in)  is  of  the  order  of  -0.3136  GPa/mm 
(-1.16  x  106  psi/in) 

Figure  3.2  shows  the  progression  of  the  ignition  point  through  the 
bed;  this  locus  is  defined  as  the  flame  front.  As  mentioned  earlier  the 
total  burn  time  is  somewhat  longer  than  seen  in  the  Krier-Kezer le  baseline, 
but  the  velocities  at  various  x  locations  are  approximately  the  same.  The 
pressure  front  as  shown  in  this  figure  was  defined  as  being  the  locus  of 
points  midway  between  the  peak  of  the  continental  divide  and  the  ambient 
pressure  level  in  front  of  the  buildup.  Using  this  definition,  the  loca¬ 
tion  of  the  pressure  front  at  60  psec  has  been  indicated  on  Fig.  3.1  with 
an  *. 

A  comparison  of  the  particle  and  gas  temperature  development  is  shown 
in  Fig.  3.3.  Again  the  continental  divide  structure  seen  in  the  gas  tem¬ 
perature  presented  in  Krier-Kezer le  is  seen  here.  The  severe  gas  tempera¬ 
ture  spike  seen  at  80  nsec  is  due  more  to  the  encounter  of  pressure  re¬ 
flection  with  the  end  wall  than  to  some  drastic  change  in  the  progression 
of  the  deflagration  wave.  It  should  be  noted  that  no  disassociation  of 
the  gases  is  assumed  in  these  calculations.  If  real  gas  properties  were 
being  assumed,  temperatures  of  this  magnitude  (i.e.,  7500°R  or  greater) 
would  certainly  cause  the  gases  to  disassociate,  reducing  the  high  values 
predicted.  The  particle  temperature  profiles  exhibit  a  structure  similar 
to  those  presented  in  Ref.  2,  i.e.,  no  continental  divide,  but  the  magni¬ 
tude  of  the  temperature  at  later  times  during  the  convective  flow  sequence 
was  brought  about  by  the  extremely  small  volume  taken  up  by  the  solid 
phase  (porosities  of  0.96  or  greater)  and  the  high  temperature  of  the  sur¬ 
rounding  gas.  Under  such  conditions,  even  the  restricted  heat  transfer 
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Fig.  3.2  Locus  of  ignition  (flame)  front  and  pressure  front,  the  latter 
derived  from  pressure  versus  distance  distribution. 
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from  the  gas  (10  percent  after  ignition)  causes  a  tremendous  increase  in 
the  temperature  of  the  particles.  Since  it  was  felt  that  this  type  of 
phenomenon  was  a  product  of  the  calculation  procedure  and  not  in  fact  real, 
subsequent  cases  remove  all  heat  transfer  at  the  higher  levels  of  porosity 
(<p  -  0.96).  Even  with  the  extremely  high  particle  temperature  exhibited  in 
this  case,  the  progress  of  the  convective  process  has  not  been  affected 
since  the  particle  phase  at  this  point  is  an  insignificant  fraction  of  the 
flow. 

The  porosity  distribution  (Fig.  3.4)  exhibits  a  contour  which  is 
characterized  by  a  zone  in  which  the  material  is  compressed  to  a  porosity 
lower  than  the  original  value.  The  point  of  lowest  porosity  is  always 
located  ahead  of  the  peak  pressure  point,  but  it  should  be  noted  that  the 
compressed  zone  does  not  prevent  the  point  of  first  ignition  from  moving 
through  this  zone.  This  ignition  point  was  always  found  to  occur  at 
relatively  low  gas  temperature,  typically  less  than  1000°R.  It  would  seem 
that  the  bed  would  need  to  be  much  more  severely  compacted  to  prevent  any 
hot  gases  from  "seeping  through"  before  a  significant  pressure  can  be  built 
up  behind  the  ignition  front.  This  seems  to  be  a  prerequisite  for  transi¬ 
tion  to  detonation. 

3.3  Variation  in  the  Boundary  Conditions  and  Grid  Spacing 

The  effects  of  the  choice  of  initial  conditions  and  boundary  condi¬ 
tions  on  a  fluid  mechanical  problem  are  known  to  be  significant  and  this 
situation  is  no  exception.  In  previous  analyses  by  Krier-Gokhale  [1]  and 
Krier-Kezerle  [2]  the  boundary  conditions  were  chosen  in  such  a  way  that 
gradients  at  the  wall  were  not  constant  with  time. 
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Fig.  3.4  Porosity  distribution  history.  The  (*)  indicates 
the  location  of  the  ignition  front. 


It  was  decided  here  that  the  problem  would  be  specified  to  assure 
that  at  every  instant  the  flow  gradients  would  be  zero  at  least  at  one  "wall 
Such  zero-gradient  boundary  conditions  are  obviously  valid  from  an  ignition 
and  flame  spreading  sequence  which  begins  in  the  middle  of  a  packed  bed. 

One  can  then  analyze  the  problem  in  one  direction  only  specifying  the 
"middle"  as  x  =  0.  Using  this  approach  effectively  removes  at  the  wall 
the  3/9x  terms  that  appear  in  the  six  conservation  equations.  By  symmetry 
arguments  the  gas  and  particle  velocities  are  zero  at  the  ignition-end  wall 
and  are  zero  at  the  far  wall  because  of  the  imposed  impenetrable  boundary. 
Due  to  the  imposed  boundary  conditions,  care  must  be  taken  that  the  initial 
conditions  (t  =  0)  satisfy  all  boundary  conditions. 

Having  established  the  conditions  for  which  zero-gradient  boundary 
conditions  could  be  used,  a  comparison  between  the  use  of  these  boundary  con 
ditions  with  the  non-zero  boundary  conditions  used  in  the  baseline  is  pre¬ 
sented  in  Fig.  3.5.  The  sensitivity  of  the  flow  to  this  change  appears  to 
be  due  to  the  fact  that  the  zero  gradient  assumption  effectively  removes 
all  influence  of  any  gradient  terms  from  the  first  and  last  grid  space. 

This  is  especially  significant  during  the  early  development  of  the  defla¬ 
gration  since  it  is  these  gradient  terms  in  the  conservation  equations 
which  give  the  gases  the  initial  impetus  to  move  forward. 

Since  the  initial  conditions  are  the  same  in  both  cases,  the  mass 
generation  rate  will  at  first  be  approximately  the  same.  But  as  the  defla¬ 
gration  process  continues,  a  slight  change  does  develop  as  shown  in  the 
results  given  in  Table  3.2.  This  information  indicates  that  when  the  igni¬ 
tion  point  reaches  a  certain  location,  the  zero-gradient  case  allows  burn¬ 
ing  for  a  longer  period  of  time  due  to  the  slower  progression  of  the  flame 
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through  bed.  Thus  the  propellant  has  had  time  to  generate  more  gas  in  a 
fixed  amount  of  space  which  accounts  for  the  higher  gas  pressure. 


TABLE  3.2 


Ignition  Point  at  50  mm  (1.95  in) 

Time 

Pg 

V 

max 

t 

(Msec) 

(Kpsi) 

(mm/psec) 

Zero  gradient  58.9 

291.5 

1.26 

Non-zero  gradient  53.9 

213.6 

1.55 

The  choice  of  the  grid  spacing  in  any  numerical  integration  scheme 
is  as  of  much  concern  as  the  choice  of  boundary  conditions.  Once  a  bed 
length,  L,  is  specified  one  would  wish  to  determine  the  minimum  number  of 
grid  spaces  which  will  adequately  handle  all  representations  of  spacial 
derivatives,  in  order  to  reduce  computation  costs  and  the  number  of  in¬ 
tegrations  required.  Table  3.3  shows  that  for  the  range  of  Ax  used  (which 
was  found  to  be  large  enough  to  observe  the  overall  trends),  there  is  not 
much  difference  in  the  final  gas  pressure  and,  to  a  certain  extent,  in  the 
total  time  needed  to  ignite  the  entire  bed.  However,  the  total  number  of 
integrations  needed  to  complete  the  calculations  is  almost  inversely  pro¬ 
portional  to  Ax,  as  indicated  in  the  table.  Figure  3.6  was  then  con¬ 
structed  in  an  attempt  to  determine  at  what  point  a  diminishing  return  was 
being  obtained  from  decreases  in  Ax.  As  can  be  seen,  there  is  not  much 


TABLF.  3.3 


Ignition  Point 

at  76.2  mm  (3.0  in) 

Ax 

t 

Pg 

^max 

Integrations 

(usee) 

(Kpsi) 

0.100 

76.76 

383.0 

84 

0.050 

72.57 

385.2 

145 

0.025 

71.29 

379.8 

280 

variation  between  a  Ax  of  0.050  and  0.025.  This  led  to  a  decision  to  us  a 
Ax  of  0.050  in  all  subsequent  calculations  since  the  results  were  close  to 
those  obtained  using  smaller  values  of  Ax.  From  to  to  time  this  comparison 
was  repeated  for  other  flow  situations  to  insure  that  the  results  were  es¬ 
sentially  insensitive  to  the  chosen  Ax  value. 

3.4  Variation  in  Propellant  Properties 

Experimental  investigations  of  the  DDT  phenomenon  in  porous  beds 
carried  thus  far  are  in  basic  agreement  concerning  the  necessity  for  a 
rapid  pressure  buildup  during  the  deflagration  phase.  Figures  3.7  and  3.8 
provide  a  comparison  of  the  pressure  buildup  and  flame  spreading  rate  for 
two  different  propellants  (both  assumed  incompressible).  The  propellant 
labeled  "energetic  composite"  has  properties  of  a  HMX  modified  propellant 
and  reflects  the  type  of  input  data  throughout  this  investigation.  On  the 
other  hand,  the  data  used  for  the  propellant  labeled  "single  base"  is  more 
typical  of  the  nitrocellulose  propellants.  As  can  be  seen,  the  energetic 
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Fig.  3.-  effects  of  propellant  energy  content  on  the 
gas  pressure  distribution  history. 
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propellant  has  nearly  burned  through  the  entire  bed  after  70  psec  having 
generated  a  peak  pressure  of  approximately  350  Kpsi.  At  this  same  instant, 
the  single  base  propellant  has  ignited  only  half  of  its  bed,  generating  a 
peak  gas  pressure  of  only  25  Kpsi.  This  indicates  the  expected  result  that 
conditions  which  might  lead  to  shock  formation  and  eventual  detonation  re¬ 
quire  high  energy  and  relatively  rapid  burning  rate  propellant  for  a  given 
bed  solids  loading  and  granulation. 

However,  it  should  be  noted  that  the  pressures  predicted  for  the 
energetic  compositie  in  Fig.  3.7  are  much  greater  than  those  reported  in 
certain  experiments  with  porous  explosives,  in  which  a  DDT  occurred.  (See 
the  comment  by  R.  D.  Bernecker  to  the  results  reported  by  Krier  and  Kecerle 
in  the  Seventeenth  Combustion  Symposium.)  The  assumption  of  an  incom¬ 
pressible  solid  is  one  reason  for  such  relatively  high  pressures.  A  pos¬ 
sible  remedy  for  this  situation  is  a  reduction  in  the  volume  taken  up  by 
the  solids  thus  leaving  a  larger  volume  into  which  the  gas  can  expand. 

This  can  be  accomplished  by  assuming  that  the  particles  are  compressible. 

As  mentioned  earlier,  a  modified  Tait  equation  gives 


[3?  ni/3 
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o  K 


The  effect  of  this  type  of  state  equation  on  the  particle  density,  and  thus 
on  its  radius,  is  shown  in  Table  3.4  for  Kq  =  2.00  x  105  psi  (a  soft  plastic) 
In  order  to  conserve  the  particle  mass'  as  the  density  increases  the  par¬ 
ticle  sice  must  decrease  according  to  the  ratio,  (r/r  )  =  (p/p^1/3. 
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TABLE  3.4 


P  (psi) 

P/P0 

r/r 

o 

SO , 000 

1.205 

0.5714 

100,000 

1..357 

0.4000 

150,000 

1.481 

0.3077 

Figure  3.9  shows  the  drastic  reduction  in  the  gas  pressure  that  is  pre¬ 
dicted  with  the  inclusion  of  this  particle  state  equation.  At  any  particu¬ 
lar  time,  the  magnitude  of  that  pressure  is  nearly  cut  in  half.  Sur¬ 
prisingly,  Fig.  .3.10  shows  that  there  is  no  significant  reduction  in  the 
time  needed  to  completely  ignite  the  bed.  This  is  caused  by  the  fact  that 
the  decrease  in  particle  radius  allows  an  increase  in  bed  porosity  which 
allows  more  hot  gases  to  move  farther  into  the  bed  in  a  shorter  period  of 
time.  However,  the  lower  gas  pressures  and  pressure  gradients  caused  by 
this  increase  in  porosity  have  exactly  the  opposite  effect  in  that  the  ac¬ 
celeration  terms  in  the  gas  momentum  equation  are  reduced  yielding  lower 
forward  gas  velocities.  The  net  effect  of  these  two  phenomena  is  to 
effectively  cancel  each  other  out,  thus  causing  little  change  in  the  flame 
front  curve. 

The  significant  lowering  of  the  gas  pressure  also  has  its  effect  on 
the  temperatures  within  the  bed.  Figure  .3.11  shows  the  reduction  of  the 
gas  temperature  to  much  more  reasonable  levels  when  compared  to  results  ob¬ 
tained  from  a  comparable  incompressible  case.  Again,  the  high  temperatures 
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Fig.  3.10  Flame  front  locus  as  a  function  of  the  solid  phase  compressibility. 
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Fig.  3.11  Gas  temperature  comparison  as  a  function  of  the 
solid  phase  compressibility. 
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at  the  far  end  wall  are  due  to  a  reflection  from  that  wall  rather  than 
from  some  significant  change  in  the  burning  process.  The  particle  tem¬ 
peratures  exhibited  in  the  compressible  case  were  also  much  lower  due  to 
the  reduction  in  the  heat  transfer  from  the  gas. 

Figures  3.12  and  3.13  allow  for  a  comparison  of  the  effects  that  two 
different  values  of  the  initial  bulk  modulus  have  on  the  properties  of  the 
system.  The  value  of  =  5.077  x  105  psi  was  obtained  from  Traver  et  al . 
[8]  who  were  investigating  the  DDT  phenomenon  in  homogeneous  solid  ma¬ 
terials.  The  second  value  of  the  bulk  modulus  used  (Kq  =  2.00  x  105  psi) 
was  obtained  from  the  sound  speed  calculation  for  the  solids  which  was 
used  by  Krier  and  Kezerle  as  part  of  their  stability  criterion.  This 
value  was  found  to  be  typical  of  a  soft  plastic  type  material. 

It  should  be  noted  in  Fig.  3.13  that  for  the  Kq  =  »  (incompressible) 
case,  the  entire  bed  did  not  ignite.  This  failure  was  due  to  the  presence 
of  numerical  instabilities  within  the  integration  scheme.  This  type  of 
problem  was  found  to  occur  in  regions  of  very  steep  gradients  which  for 
this  case  was  brought  about  by  a  slight  increase  in  the  magnitude  of  the 
interphase  drag.  This  sort  of  problem,  and  measures  which  can  be  used  to 
correct  it,  are  the  subject  of  the  next  section. 

3.5  Numerical  Smoothing  Techniques 

The  onset  of  severe  gradients  within  the  bed  often  caused  a  situa¬ 
tion  in  which  the  numerical  integration  scheme  exceeded  its  capacity  to 
provide  reasonable  answers.  The  result  was  usually  an  overshoot  of  the 
absolute  temperature  or  pressure  into  the  negative  (and  thus  physically 
impossible)  region.  Despite  the  extensive  use  of  artificial  damping 
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techniques  in  the  literature  for  similar  problems,  it  was  felt  that  the 
inclusion  of  this  type  of  damping  would  mask  important  information  concern¬ 
ing  the  onset  of  a  transition  to  detonation.  For  this  reason,  no  specific 
damping  terms  were  included  in  the  integration  scheme. 

However,  it  was  found  that  a  reduction  in  the  magnitude  of  the  in¬ 
terphase  drag  term  would  allow  progress  to  occur  for  cases  that  would 
otherwise  be  halted  due  to  the  overshoot  described  above.  This  was  due  to 
the  fact  that  less  drag  allowed  a  higher  gas  velocity  to  occur  at  a  certain 
level  of  porosity  which,  in  turn,  allowed  gases  hot  enough  to  cause  igni¬ 
tion  to  penetrate  much  more  deeply  into  the  bed.  Once  the  particles  had 
been  compacted  to  a  porosity  somewhat  below  that  which  caused  the  numerical 
oscillations  at  the  higher  drag  level,  the  hot  gases  would  again  be  unable 
to  move  into  the  bed  resulting  in  ever-steepening  gradients  and  eventual 
numerical  integration  failure. 

This  situation  is  i! lstrated  quite  clearly  in  Figs.  3.14  and  3.15. 
Figure  3.14  shows  the  ability  of  the  hot  gases  to  move  farther  and  faster 
into  the  bed  as  the  interphase  drag  is  multiplied  by  a  constant  factor,  C^. 
The  pressure  curves  (Fig.  3.15)  also  indicate  the  relaxation  of  the  bed  as 
the  drag  is  reduced.  The  results  presented  in  these  two  graphs  would  seem 
to  indicate  that  this  factor  is  acting  as  a  pseudo-damping  coefficient. 
The  constant  was  originally  used  to  allow  for  variations  in  the  gas-par¬ 
ticle  interaction  due  to  the  uncertainty  in  the  drag  relation,  which,  at 
the  high  Reynold's  numbers  encountered  in  this  type  of  flow,  was  being  ex¬ 
tended  outside  its  evaluated  range.  Until  further  experiments  can  extend 
the  range  of  validity  of  the  interphase  drag  equation,  the  use  of  a  Cp 
type  coefficient  must  remain  a  legitimate  option. 
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Flame  front  locus  as  a  function  of  bed  permeability 
coefficient,  Cn. 
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Use  of  this  factor  did  point  to  the  increased  range  of  problems 
which  could  be  investigated  if  some  sort  of  artificial  smoothing  was  im¬ 
plemented.  The  fear  of  loosing  important  data  subsided  as  it  became  more 
and  more  apparent  that  the  instabilities  seemed  to  be  resulting  from  the 
integration  scheme  itself  and  not  from  some  quality  of  the  flow.  To  this 
end,  a  damping  equation  of  the  form 
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was  utilized,  where,  the  smoothing  percent,  AVG,  must  lie  within  the  range 

0.0  >  AVG  >  0.50  (3.5) 


This  allows  up  to  SO  percent  of  the  time  derivatives  on  either  side  of  the 
term  being  worked  on  to  be  used  in  obtaining  a  smoother  result.  Figure  3.16 
shows  the  effects  of  increasing  the  amount  of  smoothing  on  the  progression 
of  the  flame  front  through  the  bed.  For  five  percent  smoothing,  the  gra¬ 
dients  remained  too  steep  to  be  handled  by  the  routine.  However,  for 
smoothing  of  roughly  10  percent  or  more,  the  integration  scheme  was  able  to 
successfully  predict  results  through  the  entire  bed. 

This  increased  success  is  not  without  its  penalties.  Figure  3.17 
shows  the  resulting  pressure  profiles  for  various  amounts  of  smoothing.  As 
can  be  seen, an  increase  in  the  amount  of  smoothing  has  the  elfect  of  smear¬ 
ing  out  any  sharp  gradients  to  the  point  where  useful  information  cannot 
always  be  obtained.  Recognizing  this  fact,  the  amount  of  smoothing  was 
kept  as  low  as  possible,  typically  less  than  10  percent. 


Fig.  3.16  Flame  front  locus  variation  due  to  numerical  smoothing 
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3 . 6  Discontinuous  Jump  Conditions  in  the  Gas  Generation  Function 

As  discussed  in  Chapter  One,  in  order  to  develop  a  strong  pressure 
shock  and  an  expected  abrupt  change  in  the  flame  velocity,  one  must  model 
a  rapid  increase  in  one  of  the  components  of  the  mass  generation  term,  T, 
at  some  assumed  critically  high  pressure  or  in  some  region  of  excess  solids 
compaction.  The  term  I'  (mass/volume-time)  for  spherical  particles  burning 
on  their  outer  surfaces  only,  is 
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(1  -  <>)Pprp 


(3.4) 


with 


-1—  =  surface  to  volume  ratio  for  spheres 

rp 
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propellant  density 


r 

P 


propellant  burning  rate 


[bPn] 


This  expression  allows  for  a  number  of  possible  changes  in  individual 
components  or  combinations  of  components  at  some  preselected  critical 
cond i t i on . 

figure  3.1s  shows  the  results  from  two  of  these  possibilities. 

Case  f .» 4  represents  the  unperturbed  prediction  of  the  flame  front  locus. 
Ca<e  (b)  illustrates  the  enhancement  of  the  flame  moving  through  the  bed 
by  simply  increasing  f  throughout  the  run.  In  this  case  the  enhancement 
was  the  result  of  allowing  the  propellant  burning  rate  to  be  augmented  by 
the  magnitude  of  the  propellant  temperature  increase.  The  third  case 
(case  ci)  made  use  of  an  increase  in  the  pressure  exponent  at  some  im¬ 
posed  critical  condition,  in  this  case  where  the  gas  pressure  exceeds 
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.18  Flame  front  locus  as  a  function  of  burning  rate  equations 
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50  Kpsi.  The  shape  of  the  flame  front  is  now  beginning  to  resemble  the 
data  of  Bernecker  and  Price  [6]  (discussed  in  Chapter  One),  although  the 
ignition  front  velocity  is  still  less  than  those  reported  in  Ref.  6 
This  case  did  not  proceed  to  completion  due  to  the  same  type  of  numerical 
instabilities  described  in  the  previous  section. 

Tn  order  to  limit  the  effect  of  the  nonlinear  source-term  increase, 
which  causes  the  numerical  integration  breakdown,  case  (c)  was  repeated 
by  assuming  only  a  five  percent  increase  in  the  exponent  n.  As  expected, 
the  integration  proceeded  normally.  As  the  pressure  history  profiles  shown 
in  Pig.  3.19  indicate,  there  is  a  significant  increase  in  the  level  of  the 
gas  pressure  shortly  after  the  propellant  burning  rate  is  increased.  The 
effect  of  this  pressure  wave  increase  on  the  motion  of  the  flame  front  is 
shown  in  Fig.  5.20.  ft  should  be  noted  here  that  while  the  gas  pressure 
exceeds  the  critical  level  at  40  Msec,  the  break  in  the  slope  of  the  flame 
front  does  not  occur  for  another  10  Msec.  However,  when  the  location  of 
the  pressure  front  is  defined  as  the  point  of  maximum  rate  of  change  of 
slope  on  the  pressure  time  plot  (Pig.  3.21),  as  was  done  by  Bernecker  and 
Price  [6],  it  can  be  seen  that  the  pressure  wave  moves  ahead  of  the  flame 
at  about  the  time  of  the  change  in  the  burning  rate  exponent.  Despite 
the  lack  of  high  flame  velocities  after  the  transition,  which  may  be  due  to 
the  relatively  low  initial  solids  loading  treated  here,  this  flame  front/ 
pressure  front  exhibits  a  remarkable  similarity  to  the  experimental  re¬ 
sults  of  Bernecker  and  Price. 

Figures  3.22  and  5.23  (flame  front  loci  plots)  show  similar  be¬ 
havior  when  changes  are  made  to  the  other  components  of  the  gas  generation 
term.  Changes  in  the  burning  rate  proportionality  constant,  b,  do  not  show 
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Fig.  3.19  Pressure  distribution  history  as  a  function 
of  burning  rate  exponent  change. 
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sure-time  variation  at  various  locations  in  bed  used  to 
in  pressure  front  shown  in  Fig.  3.20. 
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as  dramatic  a  change  in  the  slope  of  the  flame  front  simply  because  the 
assumed  increased  values  of  b  were  not  as  large  as  the  increase  due  to  the 
pressure-power  law  function.  An  assumed  abrupt  change  in  the  solid  den¬ 
sity  at  some  prescribed  high  pressure  produces  results  similar  to  those 
found  when  the  addition  of  a  state  equation  was  made,  namely,  a  reduction 
in  the  slope  of  the  flame  front  instead  of  an  increase. 

3. 7  Modification  of  the  Solid  Phase  Momentum  Equation  and 
the  Introduction  of  Stress  Wave  Motion 

The  previous  section  introduced  a  possible  mechanism  which  yields 
results  supportive  of  the  first  DDT  scenario  described  in  Section  1.2. 

But  in  all  cases  described  to  this  point  there  was  no  evidence  of  a  com¬ 
pression  wave  moving  through  the  unignited  portion  of  the  bed.  This  meant 
that  conditions  for  the  second  DDT  scenario  were  not  being  modeled.  Modi¬ 
fications  necessary  to  obtain  this  type  of  physical  effect  included  the 
introduction  of  a  new  compaction  resistance  law,  a  modification  of  the 
solid  phase  momentum  equation  and  the  assumption  that  a  solid  phase  stress 
wave  will  move  through  the  unignited  part  of  the  bed  at  the  solid  phase 
sound  speed  behind  the  wave.  The  details  of  each  of  these  new  assumptions 
are  presented  below. 

In  Chapter  Two,  problems  encountered  with  the  particle-particle 
interaction  law  developed  by  Kuo  and  Summerfield  [11]  were  discussed  and 
a  substitute  relation  based  on  the  minimum  allowable  compaction  was  out¬ 
lined.  This  new  compaction  resistance  law  was  to  be  used  in  association 
with  a  n -w  solid  phase  momentum  equation  which  took  the  final  form 


♦  [1-fWKF  -  (1-4))  ~  (f(4))F)]  (3.5) 


where,  as  defined  previously. 


p2  =  a-4))pp 


f(<t>)  = 


(3.6) 


(3.7) 


<f>UL  =  porosity  above  which  the  particles  do  (3.8) 

not  touch 


4>ll  =  minimum  allowable  porosity  based  on  (3.9) 

the  geometry  of  individual  particles 


fWF  =  Pp 


(3.10) 


Figures  3.24  and  3.25  show  the  the  effects  on  the  flame  front  locus  and 
the  pressure  profiles  of  the  addition  of  these  two  modifications  for  various 
values  of  the  exponent  l.  For  i  set  equal  to  zero,  f(<J>)  is  equal  to  unity 
for  all  values  of  porosity  less  than  <4^.  This  means  that  once  the  par¬ 
ticles  are  touching,  they  cannot  be  compacted  further  since  they  cannot  be 
accelerated  forward  and  all  the  force  exerted  by  the  drag  is  taken  up  by  a 
localized  stress  on  the  particles.  As  the  value  of  l  is  increased,  more 
and  more  of  the  drag  force  is  used  to  accelerate  the  particles  forward. 

This  results  in  the  flame  front  moving  into  the  bed  at  a  faster  rate  since 
the  particles  are  being  moved  forward  at  a  faster  rate.  This  also  means 


i 


resistance  function. 


.25  Pressure  distribution  as  a  function  of  particle-particle  resistance 
function  (at  time  of  82  psec). 
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that  there  is  additional  volume  available  for  the  gases  behind  the  flame 
front  which  results  in  a  lower  gas  pressure.  Both  of  these  events  are 
illustrated  in  Figs.  3.24  and  3.25. 

A  comparison  of  a  case  using  the  Kuo-Summerfield  formulation  for 
the  compaction  resistance  with  a  case  using  the  equation  described  above 
is  shown  in  Figs.  3.26  and  3.27.  For  this  comparison  a  value  of  unity  for 
the  exponent  H  was  chosen  since  it  represents  a  simple  linear  decrease  in 
the  acceleration  caused  by  the  drag  and  forces.  As  can  be  seen,  for 
this  case  the  Kuo-Summerfield  law  allows  a  bit  more  compaction  to  take 
place  which  in  turn  yields  a  slightly  higher  flame  spreading  rate  and  a 
lower  gas  pressure.  Thus,  if  a  value  o^  £  somewhat  greater  than  one  were 
assumed,  a  set  of  results  similar  to  those  produced  by  the  Kuo-Summerfield 
equation  would  be  obtained  without  the  questions  raised  in  Chapter  Two. 
Even  with  the  ability  to  produce  similar  results  at  these  relatively  high 
porosities,  the  results  would  become  different  as  <fi  is  reached  since  the 
Kuo-Summerfield  law  allows  the  further  compaction  of  the  particles  below 


Having  now  defined  a  particle  pressure  in  Eq.  (3.10)  as 

P  =  f(<t>)F  (3.10) 

it  is  possible  to  define  a  sound  speed  through  the  solid  phase  based  on  the 
logic  developed  in  Appendix  A.  This  information  can  now  be  used  in  con¬ 
junction  with  an  application  of  basic  continuum  mechanics  to  define  the 
velocity  with  which  a  stress  wave  will  move  through  the  unignited  bed. 

Figure  3.2S  shows  clearly  the  progression  of  this  stress  wave  through  the 
bed  far  in  advance  of  the  gas  pressure  front.  Since  it  was  assumed  that 
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this  stress  wave  could  affect  the  particle  density  in  the  same  manner  as 
the  gas  pressure,  the  particles  are  compressed  and  moved  slightly  forward 
by  this  Pp  force.  The  net  effect  is  to  increase  the  porosity  in  front  of 
the  gas  pressure  wave  allowing  hot  gases  to  penetrate  forward  much  more 
quickly  and  ignite  a  greater  portion  of  the  bed  in  the  same  amount  of  time. 
This  increase  in  the  porosity  ami  the  associated  increase  in  the  flame 
speed  are  shown  in  Figs.  3.29  and  3.30. 


However,  it  is  also  well  known  from  fundamental  continuum  mechanics 


that  the  formation  of  an  insipient  shock,  which  was  the  goal  behind  this 
particular  investigation,  will  only  result  if  the  source  of  the  compression 
waves  is  accelerating  forward.  The  case  which  is  represented  in  Figs.  3.28 
through  3.30  does  not  exhibit  this  necessary  behavior  until  the  final  stages 
of  the  burn  with  the  result  that  the  particle  stress  is  increased  in  a 
gradual  manner  until  the  bed  is  completely  ignited  with  no  sign  of  shock 
format  ion. 


A  possible  remcdv  to  this  situation  is  to  increase  the  bu"  modulus 
of  solid  phase  which  would  have  the  dual  effect  of  increasing  the  forward 
’■eloeity  of  tie  compression  waves  and  decrease  the  amount  of  compression  in 
the  solid  phase.  Increasing  the  sound  speed,  i.e.,  the  compression  wave 
speed,  would  increase  the  potential  for  shock  formation  in  a  bed  of  this 
short  length.  A  decrease  in  the  amount  of  compression  will  under  these  cir¬ 
cumstances  have  the  effect  of  slowing  down  the  progression  of  the  flatne 
front.  This  is  so  since  the  porosity  will  not  be  increased  as  much.  Cal¬ 
culations  with  an  assumed  increased  bulk  modulus  showed  that  indeed  the 
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Fig.  3.29  Porosity  distribution  comparison  for  localized 
and  non-local i zed  solid  phase  stress. 


Fig.  3.30  Flame  front  locus  comparison  for  localized  and  non-localized  solid 
phase  stress. 
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wave  source  was  not  proper  to  develop  an  insipient  shock.  The  increase  in 
the  solid  phase  bulk  modulus  does  have  the  effect  of  compressing  a  suffi¬ 
cient  number  of  particles  against  the  back  wall  to  cause  the  gas  to  be 
compressively  heated  due  to  the  decrease  in  the  volume  available  for  it. 
This  eventually  resulted  in  the  gas  obtaining  a  high  enough  temperature  to 
ignite  the  particles  before  the  arrival  of  the  deflagration  front.  How¬ 
ever,  this  was  a  result  of  the  fact  that  the  wall  allowed  no  further  motion 
of  the  particles  and  not  from  the  desired  shock  initiation. 

Based  on  the  investigations  presented  in  this  section,  it  appears 
that  there  is  a  great  deal  of  potential  in  exploring  this  latter  scenario. 
It  also  shows  that  a  great  deal  more  must  be  known  about  the  mechanics  of 
a  bed  of  closely  packed,  individual  solid  particles.  A  large  number  of 
assumptions  had  to  be  made  here  that  were  not  entirely  satisfactory  due  to 
this  lack  of  information. 
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CHAPTER  FOUR 
CONCLUSIONS 


4. 1  Summary  of  Investigations 

The  work  presented  in  this  report  represents  a  continuation  of  the 
efforts  begun  by  Kr ier-Gokhale  [1]  and  Krier-Kezer le  [2]  to  correctly  model 
the  deflagration  to  detonation  transition  phenomenon.  In  these  previous 
studies,  the  investigation  was  focused  on  determining  a  set  of  governing 
equations  which  properly  described  the  peculiar  nature  of  the  flow  and 
with  properly  incorporating  these  equations  into  a  numerical  integration 
scheme.  This  study  essentially  began  with  these  developed  codes  and  to  de¬ 
termine  the  sensitivity  of  the  various  assumptions  made  previously  in  an 
attempt  to  better  understand  the  properties  of  the  flow  which  are  necessary 
to  correctly  model  a  smooth  transition  to  detonation. 

To  this  end,  the  boundary  conditions  and  grid  spacing  were  examined 
and,  in  the  case  of  the  boundary  condition,  were  modified  to  reflect  the 
assumption  that  the  gradients  at  the  boundaries  should  be  zero.  The  di¬ 
mensions  of  the  grid  spacing  was  varied  until  a  size  was  found  which  gave 
proper  accuracy  and  minimum  computing  costs. 

Concerning  basic  properties  of  the  propellant  itself,  it  was  shown 
that  the  energetic  composite  class  of  materials  was  much  more  sensitive  to 
the  rapid  increases  in  pressure  and  in  flame  velocity  than  were  the  single 
base  propellants,  as  expected.  It  was  also  shown  that  introducing  a  solid 
phase  state  equation  to  allow  the  particles  to  compress  dramatically  im¬ 
proves  the  characteristics  of  the  flow  properties  when  compared  with  the 
results  obtained  in  experiments.  The  introduction  of  compressibility 
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allowed  for  a  significant  decrease  in  the  gas  pressure  and  temperature  to 
levels  more  in  line  with  the  values  found  in  these  experiments,  while  re¬ 
taining  approximately  the  same  flame  front  velocity  found  using  the  incom¬ 
pressible  assumption. 

In  an  attempt  to  increase  the  range  of  problem  that  could  be  ac¬ 
commodated  by  this  code,  a  smoothing  function  was  introduced  to  suppress 
numerical  oscillations  which  occurred  with  certain  input  conditions.  These 
oscillations  were  also  found  to  be  suppressed  by  a  reduction  in  the  grid 
which,  through  the  Courant-Friedrichs-Levy  stability  criterion  [12],  also 
reduced  the  time  step.  But  a  reduction  of  this  type  of  course  increased 
the  computing  time,  costs,  and  even  accuracy.  It  was  found  that  numerical 
stability  could  be  improved  by  gradually  reducing  the  time  step  below  the 
Courant-Friedrichs-Levy  criterion  as  the  gradients  became  more  severe  with¬ 
out  changing  the  grid  site.  While  this  also  increased  the  computer  time 
needed,  the  increase  was  not  significant.  This  was  due  to  the  fact  that 
although  the  number  of  integration  cycles  was  increased,  the  number  of 
steps  per  cycle  remained  the  same. 

At  this  point,  an  investigation  of  two  different  theories  concerning 
the  DDT  mechanism  was  carried  out  to  determine  what  type  of  modifications 
would  be  necessary  to  model  these  two  theories.  In  order  to  reproduce  the 
first  mechanism,  which  has  its  origins  in  experimental  work  by  Burnecker 
and  Price  [6] ,  a  modification  which  yielded  positive  results  was  to  abruptl 
increase  the  magnitude  of  the  gas  generation  term.  The  most  dramatic 
effects  of  such  a  change  were  produced  by  changing  the  exponent  in  the 
pressure  sensitive  burning  rate  law  at  some  imposed  critical  pressure  level 

The  second  mechanism,  which  had  its  origin  in  the  work  of  Ma?ek  [7], 
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was  modeled  here  by  modifying  the  solid  phase  momentum  equation  which  in 
turn  allowed  a  ’.'article  stress  to  be  defined.  This  stress  term,  coupled 
with  an  expression  for  the  solid  phase  sound  speed,  created  a  situation  in 
which  compression  waves  could  move  forward  into  the  bed  with  the  potential 
for  the  formation  of  an  insipient  shock.  While  no  such  strong  shock  for¬ 
mation  has  been  calculated  as  yet,  work  is  continuing  to  expand  and  refine 
this  approach. 

4 . 2  future  Areas  of  Investigation 

The  work  presented  in  this  report  has,  among  other  things,  served 
the  purpose  of  focusing  attention  on  various  aspects  of  the  OUT  phenomenon 
which  will  require  a  much  more  thorough  understanding  before  any  significant 
confidence  in  the  results  can  be  obtained. 

for  example,  it  is  unknown  at  this  time  if  the  fundamental  relation¬ 
ships  for  the  drag  and  heat  transfer  are  valid  at  the  high  Reynolds  numbers 
encountered  in  this  type  of  flow.  Results  presented  in  Section  3.5  show 
what  kind  of  effect  a  reduction  in  the  drag  relation  can  have  on  the  pro¬ 
gress  of  the  flow.  It  is  also  uncertain  what  effect  the  large  gas  pres¬ 
sures  have  on  the  burning  rate  of  these  types  of  propellants.  That  is, 
dynamic  burning  rate  phenomena  may  negate  the  use  of  the  rate  equation  used 
in  Chapter  Three. 

The  large  pressure  buildups  in  the  short  time  periods  which  occur  in 
this  type  of  problem  would  be  expected  to  create  some  sort  of  dynamic  com¬ 
paction  in  the  porous  bed.  But  again,  the  nature  of  this  compaction  is  not 
fully  understood  and  therefore  cannot  be  properly  modeled  in  this  code, 
further  experimentation  will  need  to  be  conducted  if  questions  of  this  type 
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are  to  be  answered. 
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Further  improvements  in  the  code  itself  can  be  obtained  with  the 
introduction  of  concepts  commonly  used  in  detonation  physics.  This  would 
include  the  use  of  a  revised  gas  generation  term  which  reflects  the  in¬ 
creased  mass  transformation  rates  which  occur  in  a  detonation  [14].  Use 
would  also  need  to  be  made  of  detonation  chemistry  which  allows  for  the 
shock  initiation  of  the  types  of  materials  being  considered  here.  An 
example  of  the  possible  use  of  this  concept  is  to  allow  an  abrupt  change  in 
the  burning  rate  exponent  as  was  done  in  Section  3.6.  However,  this  in¬ 
creased  burning  rate  probably  should  be  limited  to  what  could  be  called  the 
shock  front,  that  is,  the  exponent  increase  should  be  confined  to  that  re¬ 
gion  where  the  gas  pressure  is  increasing  at  some  critical  rate  with 
respect  to  time.  As  has  been  shown  throughout  Chapter  Three,  this  type  of 
zone  is  always  found  to  be  closely  coupled  to  the  ignition  front. 

An  indication  of  the  results  of  this  type  of  localized  burning  rate 
increase  are  presented  in  Figs.  4.1  and  4.2.  Specifically,  the  exponent  n 
is  increased  from  nQ  to  nQ  +  An,  only  when  dP/dt  >  1.5  x  10 3  psi/  sec  and 
only  if  the  pressure  gradient,  dP/dx  is  negative.  Although  these  results 
resemble  those  shown  in  Section  3.6,  the  logic  used  here  models  a  sequence 
of  events  which  is  much  closer  to  the  circumstances  known  to  occur  in  a 
detonation.  The  flame  front  velocities  which  occur  after  the  break  in  its 
slope  were  found  to  be  constant,  which  is  another  important  characteristic 
of  a  detonation.  And  finally,  the  crossing  of  the  flame  front  by  the  pres¬ 
sure  front  before  the  break  in  the  flame  front  slope  was  found  to  occur  in 
all  the  cases  presented  in  Fig.  4.1.  But  due  to  the  means  by  which  the 
reaction  zone  was  defined  (in  this  case,  the  restriction  that  dP/dx  <  0), 
the  pressure  front  eventually  dropped  below  the  flame  front  for  the  higher 
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Flame  front  locus  resulting  from  an  abrupt  change  (at  the  critical 
pressure  gradient)  in  the  burning  rate  exponent.  (Limited  to  the 
localized  region  of  the  pressure  front.) 
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Fig.  4.2  Pressure  distribution  history  as  a  function  of 
burning  rate  exponent  change.  (Limited  to  the 
localized  region  of  the  pressure  front.) 
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values  of  An  as  the  burn  progressed.  This  was  caused  by  the  fact  that  in¬ 
creases  in  An  resulted  in  a  steepening  of  the  pressure-distance  curve  which 
in  turn  caused  the  reaction  cone  to  become  narrower  in  extent  and  thus  less 
powerful  in  its  motion  through  the  bed.  This  problem  can  be  rectified  by 
using  a  somewhat  more  complicated  means  of  defining  the  reaction  zone  which 
will  allow  the  thickness  of  the  zone  to  be  maintained  at  a  desired  level 
independent  of  the  pressure-distance  curve. 

It  should  be  noted  that  in  all  the  cases  presented  in  this  report, 
the  initial  solids  loading  is  only  60  percent  (porosity,  <p ,  =  0.40)  while 
tiie  experiments  of  Burtiecker  and  Price  [6]  as  well  as  those  of  Calzia  and 
Carabin  [9]  had  initial  solids  loading  ranging  approximately  from  80  to  90 
percent.  Appendix  B  shows  that  an  initial  solids  loading  of  74  percent  is 
the  maximum  attainable  with  uniform  spheres,  which  is  a  basic  assumption  in 
this  code.  To  treat  higher  values  of  solids  loading  one  must  assume  that  a 
bed  of  multisized  particles  is  being  used.  This  leads  to  a  number  of  dif¬ 
ficulties  concerning  the  proper  method  of  partitioning  flow  properties  such 
as  the  interphase  drag  and  heat  transfer.  (For  a  discussion  of  these  dif¬ 
ficulties  and  their  potential  solution  see  the  Ph.D.  thesis  of  S.  S.  Gokhale 
[Id].) 

Another  aspect  which  should  be  looked  into  both  as  a  means  of  im¬ 
proving  this  code  in  particular  and  for  solving  this  DDT  problem  in  general 
is  to  review  again  the  numerical  scheme  itself.  The  full  effects  of  the 
numerical  smoothing  function  i ntroduced  in  Section  3.S  have  not  been  ex- 
hustively  studied  and  refinements  are  certainly  not  out  of  the  question. 

The  use  at'  a  "rezoning"  technique  in  which  either  the  grid  spacing  or  the 
time  step  or  both  are  decreased  in  regions  of  severe  gradients  may  prove 
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to  be  very  valuable  in  obtaining  better  accuracy  for  the  same  or  even  lower 
computing  time  arul  costs.  Reasonable  solutions  may  not  be  acquired  from 
any  of  these  suggestions  until  larger  computers  become  available  since 
■•’.any  of  the  calculations  encountered  thus  far  entailed  calculating  small 
differences  of  large  numbers. 

Finally,  it  should  again  be  pointed  out  that,  although  a  number  of 
difficulties  and  constraints  have  been  discussed  in  this  section,  many  of 
the  problems  would  not  be  recognized  as  such  without  the  progress  that  has 
been  made  in  this  work.  Further  research  into  these  problems  will  certainly 
make  progress  towards  a  solution  to  the  overall  prediction  of  deflagration 
to  detonation  transition. 
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APPENDIX  A 


PARTICLE  BULK  MODULUS  AND  SOUND  SPEED 


The  general  equation  for  the  sound  speed,  c,  in  a  solid  is  usually 
defined  as 


c 


(A.  1) 


where  K  is  the  bulk  modulus.  But,  assuming  nothing  is  known  about  the 
variation  of  K,  this  equation  would  appear  to  force  the  sound  speed  to 
decrease  as  the  particle  density  increases  which  seems  inconsistent  with 
the  normal  occurrence  of  sound  speed  varying  in  the  same  manner  as  the 
density. 

This  indicates  that  the  bulk  modulus  also  varies  with  density  which 
can  be  determined  from  the  definition  of  the  bulk  modulus  once  a  state 
equation  has  been  specified.  It  has  been  assumed  in  this  work  that  the 
particles  obey  a  modified  form  of  the  Tait  equation 
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Since  the  bulk  modulus  is  defined  as 
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the  Tait  equation  is  the  only  expression  necessary  to  determine  its  varia¬ 
tion  with  pressure  or  density.  First,  it  should  be  noted  that  for  a 
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fixed  mass 

P-P  P-P 


o  o 


Using  this,  the  bulk  modulus  can  be  written  as 


Substituting  this  into  (A. 3b) 
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This  can  now  be  put  back  into  (A.l)  to  find  the  actual  variation  in  sound 
speed  with  respect  to  particle  density: 


c 


1/2 


(A.  6) 


This  equation  now  shows  sound  speed  increasing  as  density  increases,  as 
would  be  expected. 
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APPENDIX  B 

POROSITY  LIMITS  FOR  SPHERICAL  PARTICLES 

The  following  calculations  are  provided  so  that  a  geometrical 
association  can  be  made  with  certain  values  of  the  porosity,  <t>.  It 
should  be  noted  that  for  uniform  spherical  particles,  case  (A)  represents 
the  lower  bound  for  <f>  while  case  (C)  is  the  upper  bound. 

Case  (A)  Face  Centered  Cubic  (See  Fig.  B.l.a) 

VT"VP  b 3-4 r  (4/3)irr3] 

*  =  ^ - P - 

/2 

For  this  configuration  r  =  -j-  b.  Thus 

4,  =  b3-(16/3)  (tt)  (/2/4)  3b3 

=  0.2595 

Case  (B)  Body  Centered  Cubic  (See  Fig.  B.l.b) 

,  _  VT'VP  b3-2(4/3)  (iT)r3 

* =  -v; - b^ 

/s 

For  this  configuration,  r  =  —  b.  Thus 
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Fig.  B.l  Relative  geometry  of  spherical  particles  at  various  porosity  levels 


Case  (Q  Simple  Monoclinic  (See  Fig.  B.l.c) 

*  -  VT'VP  -  b3-  (4/3)  (TT)r3 
V  vT  b3 

For  this  configuration,  r  =  ^ 


b3-(4/3)7rb3/8 

- F - 


=  0.4764 


Case  (D)  No  Particle  Contact  (See  Fig.  B.l.d) 


VT~VP  b3-f4/31iTr 3 
■  VT  b3 


In  this  configuration,  r 
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a+2 


-  b3-  (4/3)  (tt)  (1/ fa+21  1 3b3 

^3 


=  1  -  (4tt/3)  (l/[a+2])  3 


A  plot  of  4)  vs  a  is  provided  in  Fig.  B.2.  As  can  be  seen,  as  a  approaches 
zero  <(>  approaches  the  value  found  in  Case  (C)  as  would  be  expected. 


es  for  a  fluidized  bed 


APPENDIX  C 


CONSTITUTIVE  RELATIONS 


This  section  describes  in  more  complete  detail  the  constitutive 
relations  and  criteria  outlined  in  Chapter  Two. 

Due  to  the  selection  of  separate  energy  equations  for  each  phase  of 
the  flow  it  becomes  impractical  to  use  the  method  of  Kuo  et^  al_.  [15]  to 
define  an  ignition  temperature  based  on  the  solution  of  the  heat  con¬ 
duction  equation  for  the  particles.  Instead,  the  method  proposed  by  Krier 
et  al. [16]  in  which  a  "bulk  particle  temperature",  defined  as  the  average 
temperature  of  the  solids,  was  used  as  the  critical  variable  for  determin¬ 
ing  ignition.  Under  this  approach,  ignition  occurs  once  the  bulk  tem¬ 
perature  reaches  some  critical  TjGN  which  is  less  than  the  critical  sur¬ 
face  temperature  used  by  Kuo  et^  al_.  [15]. 

Convective  heat  transfer  between  the  hot  gases  and  the  particles 
made  use  of  Denton's  [17]  heat  transfer  coefficient 


Pg 


0.65 


K 
_2_ 
2r 


Pg|Ug-Up|W(2rt>) 


g 


0 . 7 


[Pr] 


0.33 


(C.l) 


where  K  is  the  thermal  conductivity  of  the  gas,  u  is  the  viscosity  of 
g  g 

the  gas,  and  Pr  is  the  Prandtl  number  of  the  gas  phase. 

Due  to  questions  concerning  the  alteration  of  the  boundary  layer 
around  a  burning  particle  which  in  turn  alters  the  convective  heat  trans¬ 
fer,  it  has  been  assumed  that  the  above  value  for  the  heat  transfer  co¬ 
efficient  is  reduced,  once  the  particle  is  ignited,  to  one-tenth  of  the 
value  used  if  the  particle  were  not  bunting. 

The  drag  relation  used  in  this  investigation  is  that  devised  by 
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APPENDIX  D 

ANALYSIS  OF  THE  INTERNAL  ENERGY  RELATION  WHEN  UTILIZING 
A  NONIDEAL  STATE  EQUATION 

From  simple  thermodynamics  one  can  write  the  internal  energy  as 


de  = 


fltr)  dT  +  f|S.] 
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dv 


(D.l) 


where 


e  =  the  internal  energy 
T  =  temperature 

v  =  the  specific  volume  (=  1/p) 


It  can  also  be  shown  that 
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These  last  two  equations  are  derived  in  detail  in  Lee  and  Sears  [19].  How¬ 
ever,  by  definition 


8  =  coefficient  of  volume  expansion  = 
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K  =  isothermal  compressibility  =  - 
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These  two  partial  derivatives  can  now  be  evaluated  using  the  equation  of 
state  since  they  contain  only  state  variables. 

This  analysis  will  proceed  without  regard  for  the  exact  form  of  B 
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in  order  to  keep  the  results  as  general  as  possible.  In  order  to  analyze 
Eq.  (D.4),  the  state  equation 


p  R  T 

p  =  S  §  g 

g  1-P  B 
g  v 


(D.6) 


is  differentiated  on  both  sides  with  respect  to  T  holding  P  constant 
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yielding 
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Solving  this  equation  for 


3v 


3T 


snd  dividing  both  sides  by  v  will 
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result  in  the  necessary  form  for  6 
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To  solve  for  k,  both  sides  of  Eq.  (D.6)  are  first  differentiated  with 
respect  to  while  holding  constant,  resulting  in 
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Solving  this  equation  for 
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and  dividing  both  sides  by  -v  yields 
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Now,  substituting  the  results  from  Eqs.  (D.8)  and  (D.10)  into 
Eq.  (D.3)  yields 
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Thus,  regardless  of  the  exact  form  of  Bv,  Eq.  (D.3)  will  always 
be  zero  and 

de  =  cvdT 

only. 

c 

However,  the  ratio  of  specific  heats,  y  =  — ,  will  not  be  a 
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constant  for  the  gas  that  obeys  Eq.  (D.6).  Again,  utilizing  Lee  and 
Sears  [19]  one  obtains 
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